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Abstract 

We study the non-canonical symplectic structure, or iC-symplectic structure inherited by the 
charged particle dynamics. Based on the splitting technique, we construct non-canonical symplectic 
methods which is explicit and stable for the long-term simulation. The key point of splitting is 
to decompose the Hamiltonian as four parts, so that the resulting four subsystems have the same 
structure and can be solved exactly. This guarantees the iC-symplectic preservation of the numerical 
methods constructed by composing the exact solutions of the subsystems. The error convergency of 
numerical solutions is analyzed by means of the Darboux transformation. The numerical experiment 
display the long-term stability and efficiency for these methods. 
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I. INTRODUCTION 


In the application of magnetized plasmas, the motion of charged particles under the influ¬ 
ence of electromagnetic fields is the most fundamental process. The long-term simulation of 
the dynamics with better accuracy is signihcant for various multi-timescale problems, such 
as particles in tokamaks. It is a recent successful development to apply the idea of geometric 
numerical methods in simulating the dynamics of charged particle equations. In QQ, 
according to the volume-preserving nature of the charged particle system the (high-order) 
numerical methods have been constructed in order to preserve the volume in phase space. 
Compared with the standard methods such as the fourth order Runge-Kutta method, this 
kind of volume-preserving numerical methods can bound the errors of conservative laws of 
the given system, and yield accurate numerical result over long time. 

In a given electromagnetic field E(x) and B(x), the motion of a single particle admits 
the Lorentz force law with which the dynamical system can be described as 


X = V, mv = q (E(x) -|- v x B(x)). (1) 

It is well known that the system has a canonical structure in coordinates (x, p)Q. However, 
generally classical canonical symplectic methods for this system can not be implemented 
explicitly. Notice that Eqs. ([T]) itself is separable and equipped with a non-canonical sym¬ 
plectic structure. It motivates us to construct a new kind of explicit numerical methods 
based on this property. 

With this purpose, we rewrite ([T]) as 


( —gB(x) —ml 
ml 0 



gV(p(x) 


mv 


( 2 ) 


/ 0 -S3(x) S2(x) \ 

where B(x) = S 3 (x) o -Bi(x) is defined by the magnetic field B. It is in a more 

V-S 2 (x) Si(x) 0 J 
compact form as 

i = K-\z)VH{z), (3) 


with z = [x"'", v"’"]"'", K{z) the skew-symmetric matrix on the left side of Eq. ([T]), and 
H(x, v) = mv ■ v/2 -|- g(p(x) the energy. Denote a two-form in phase space (x, v) with the 
entries of K{z) as 

D = g/2dx A B(x)dx -I- mdx A dv. (4) 
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It is clear from (jlj) that dhl = 0 for V ■ B(x) = 0. Thus, the skew-symmetric matrix K{z) 
dehnes a non-canonical symplectic structure in (x, v). This implies that the phase flow 
(pt : z I— )■ (pt(z) of system (j3]) preserves the it'-symplectic structure, which is 

^d<pt(z) A K{(pt{z))d(pt{z) = 0. (5) 

The purpose of the current work is to present a class of explicit i^-symplectic structure¬ 
preserving methods. They are constructed using the Hamiltonian splitting which decomposes 
the given system as several solvable subsystems. We organize the paper as follows. In 
section 2, we present the construction of i^-symplectic numerical methods based on splitting 
technique. We also provide the theoretical results on error of numerical solutions by Darboux 
transformation. Experiments are shown in section 3. Section 4 concludes this paper. For 
simplicity, in the following discussions the equations are all normalized. 


II. hT-SYMPLECTIC STRUCTURE-PRESERVING METHODS 

We start this section by introducing a dehnition of i^-symplectic structure-preserving 
methods. 

Definition II. 1. A numerical method z„+i = iphi'^n) applied to system (I3l) is called K- 
symplectic structure-preserving if it satisifes 

d^ph{zn) A K{^ph{zn))d^ph{z,,) = dz„ A A'(z„)dz„, (6) 

or, alternatively 

(^) 

It is clear that (E]) is a discretization of ([5]). The questions arises as to how the K- 
symplectic-preserving numerical methods can be constructed. Already it is known that 
the traditional numerical methods can not, in general, preserve the A'-symplectic structure. 
Nevertheless, possible construction can be accomplished if the concerned system can be split 
as several solvable subsystems which share the same structure as the original system. Due 
to this, we rewrite the Hamiltonian as 

.. 3 3 

i=l i=l 
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The Lorenz force system is, therefore, equivalent to 

3 

z = + (8) 

i=l 

This leads to four subsystems possessing the same structure ([5]). It is observed that all the 
subsystems can be solved explicitly. The solutions are displayed as below. 

Firstly we consider the dynamics generated by it is 


X = 0, 

V = E(x). 


( 9 ) 


It is obviously an integrable system. From the point z(t) = (x(t),v(t)), the exact solution 
flow of this subsystem denoted by z(t + h) = 0^‘^(z(t)) is 

f x(t + h) = x(t), 

( 10 ) 

y v(t + h) = v(t) + hE(x(t)). 

For a given i {i = 1, 2, 3), the subsystem associated with Hy. is 


X = ViBi, 
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( 11 ) 


i=i 

where e* is the three-dimensional vector with the i-th element being 1, and Bij is the entry in 
the ith row and jth column of B. As Vi is constant along time, the subsystem is also exactly 
solvable. From the given point (x(t),v(t)), the solution of the subsystem (ITT]) denoted by 
z{t + h) = 0f“‘(z(t)) is 


{ x(t + h) 
v(t h) 


= x(t) hvi(t)ei, 


= V 


(t) + ^ Bi 

1=1 


r‘Xi{t)+hvi (t) 


Xi{t) 


Bii{yi{t))dxi, 


More specihcally, when i = 1 the solution is 


( 12 ) 


x(t + h) = x(t) -f hvi(t)Bi, 

v(t + h) = v(t) Ff ^(x(f), v(t), h)B2 + F3^^^(x(t), v(t), h)e3, 

where B 3 {^,X 2 ,X 3 )d^ and F 3 ^^^(x,v, h) = B 2 {^, X 2 , X 3 )d^. 

We note that the integrals in Eq. flT^ should be evaluated exactly, which is feasible when 
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the magnetic field is given in terms of familiar functions. When the integrals cannot be 
evaluated exactly, we can chose a discretization of the magnetic field such that the integrals 
can be exactly evaluated under this discretization. For example, we can discretize the 
magnetic field on a Cartesian grid using piece-wise polynomials. Then the non-canonical 
symplectic structure associated with a given discretized magnetic field is preserved exactly 
by In practice, the magnetic field is often specified discretely. 

By composing the solutions in Eqs. flTOj) and (IT^ of the four subsystems, we obtain a 
/f-symplectic method for the Lorenz force system ([3]) with the accuracy of order 1 as 

o O 3 O (13) 


The iC-symplectic method up to second order accuracy can be gained by a symmetric com¬ 
position, which is 










(14) 


By definition III. 11 it is clear that the methods (lT3l) and (ITTl) are iC-symplectic because of 
the group property. Furthermore, it is possible to increase the accuracy of the numerical 
methods by various compositions js, 5|. It is obvious that the numerical methods (fT3|) and 
(lT4l) presented here are all explicit and easy to compute. 

Taking the determinate of ([7]) on both sides, we get 


det 


d^h 

dz„ 


det{K{Lph{zn))) = det(iF(z„)). 


The above equality provides det = 1 via using the definition of K{z). This implies 

that for the Lorentz force system the iF-symplectic method is also volume-preserving. 

By Darboux’s theorem, there exists a coordinate transformation z = T(y) such that a 
non-canonical system ([3]) can be turned into a canonical system 


y = J-'VF(y), 


(15) 


where H{y) = H{T{y)). Under the coordinate transformation, a iF-symplectic numerical 
method d?]) is also transformed to a symplectic method. As H{y) = H{T{y)) = H{z), under 
the conditions of Theorem IV.8.1 in we have the numerical energy error estimate for the 
A-symplectic methods of order p 

H{zn) = H{zo) + 0{hP), nh < (Ig) 
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Assume that the transformation function T is bounded, under the condition of Theo¬ 
rem X.3.1 we have the following error estimate for the numerical solution 

||z„ — z(t)|| < Cth^, for t = nh < h~^. (17) 

These results guarantee that the A'-symplectic methods can bound the energy error over 
exponentially long time. Moreover, the global solution error is restricted to linear growth, 
which is much smaller than the quadratic or exponential error growth of standard methods. 

III. NUMERICAL EXPERIMENTS 

In this section, we test the non-canonical symplectic methods presented in the above 
section by simulating the charged particle dynamics in a given electromagnetic field. We 
consider the two dimensional dynamics of a charged particle under the symmetric static 
electromagnetic field, 

10 "^ 

B = E = -^(xiei 0 : 262 ) (18) 

where R = ^Jx\ + x^- In the implementation of the numerical method (IT^ or flTTl) . it is 
noticed from ffT 2 |) that we need to calculate an integral of the magnetic field B. It usually 
can be computed explicitly. In this example, integrating B with respect to Xi gives 

pxi 

/ 53 (C X2, X^)di = f{xi, X2) - fixi - tVi,X2) 

J Xl—tvi 

with /(xi, X 2 ) = Xx^Jxx^ -t- X 2^/2 X 2 hi (xi + / 2 . 

Firstly, we illustrate the long-term behavior of a second order non-canonical symplectic 
method, compared with the traditional 4-th order Runge-Kutta method. The initial position 
and velocity are taken as xq = —62 and vq = 0.2ei -|- 0 . 1 e 2 . With the step size h = tt/IO, 
we perform the two methods over the time [ 0 , 2 x lO^h]. 

In Fig. [T]( a) and (b), numerical orbits of the two methods from the 190000-th step are 
plotted. Theoretical analysis shows that the particle’s orbit is a spiraling circle with constant 
radius, where the large circle corresponds to the drift of the guiding center, and the small 
circle is the fast-scale gyro motion. From Fig. [T](a) it is observed that after long computation 
time the 2nd order K-symplectic method still gives a correct orbit, while in Fig. [I](b) the 
4-th order RK method fails because of the dissipated numerical gyromotion. This can be 
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Figure 1. (a) Numerical orbits of the 2nd order K-symplectic method from the 190000-th step, (b) 
Numerical orbits of the 4th order RK method from the 190000-th step, (c) Energy error of the 
two methods as a function of the computing steps till the 20000-th step, (d) Solution error as a 
function of the computing steps, the step size in this sub-figure is h = 0.237r. 


explained partly by their performances at the energy. The obvions energy damping by the 
4nd order Rnnge-Kutta method is demonstrated in Fig. [^c). On the contrary, the second 
order K symplectic nnmerical method can bonnd the energy np to error (see Figs|T](c), 
|2](b)) which agrees with the theoretical error analysis. In Fig|T](d), we show the error growth 
of the nnmerical solution along with the simulation time. The step size is chosen to be 
h = 0.237r. It can be seen that the error caused by the Runge-kutta method is smaller at 
the first few steps, but grows larger than that of the iF-symplectic method rapidly. 

Next, we study the accuracy and computing amount of the iF-symplectic methods. 
Figj2](a) displays the solution error as a function of the computing cost. When the require¬ 
ment on solution error is smaller than 1%, the second order method needs less computing 
time and is more efficient than the hrst order method. 
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(a) (b) 

Figure 2. (a) Solution error as a function of the computing time, (b) Error of the numericl energy 
as a function of the step size h. 


IV. CONCLUSION 


In this paper, we have constructed numerical methods for the Lorentz force system focus¬ 
ing on its iC-symplectic structure. As there exists a Darboux transformation such that the 
A-symplectic structure can be turned into the standard symplectic structure, the error anal¬ 
ysis for the canonical numerical methods can be generalized to the A'-symplectic numerical 
methods. However, the Darboux transformation usually is not able to be expressed explic¬ 
itly. Thus, it is generally not easy to design the A'-symplectic numerical methods. For our 
study, we decompose successfully the Lorentz force system as four subsystems which can be 
solved explicitly. The resulting numerical methods by composing the exact solutions of the 
four subsystems naturally preserve the iF-symplectic structure. The numerical experiment 
shows the superior long-term performances of the new derived numerical methods. More 
theoretical and numerical analysis of the K-symplectic methods will be reported in future 
publications. The splitting idea has been generalized to construct the numerical methods 
for the Vlasov-Maxwell equations in j^. 
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